Niniejszy raport zawiera analizę danych pacjentów chorych na koronawirusa przyjętych do szpitala Tongji w Wuhan na przełomie stycznia i lutego 2020 roku. W raporcie uwględnione zostały analizy wraz z wykresami przedstawiające zależności pomiędzy różnymi cechami pacjentów a śmiertelnością. Analizowanymi atrybutami były między innymi: dane demograficzne pacjentów (płeć, wiek, etc.) oraz dane medyczne pacjentów powstałe w wyniku wykonanych na nich badaniach krwi. W raporcie przedstawione zostały również korelacje pomiędzy różnymi biomarkerami oraz dołączone zostały wykresy zmienności najbardziej kluczowych atrybutów na przestrzeni kolejnych badań. W ostatnim rozdziale raportu stworzony został klasyfikator dokonujący predykcji śmiertelności pacjenta zakażonego koronawirusem. Dokładność opisanego klasyfikatora osiągnęła poziom 95.40%.
library(readxl)
library(httr)
library(knitr)
library(plyr)
library(gridExtra)
library(tidyverse)
library(janitor)
library(scales)
library(kableExtra)
library(qwraps2)
library(rmarkdown)
library(caret)
library(reshape2)
library(plotly)
library(naniar)
library(gganimate)
Surowe dane zostały pobrane z adresu url (http://www.cs.put.poznan.pl/dbrzezinski/teaching/zed/wuhan_blood_sample_data_Jan_Feb_2020.xlsx) bezpośrednio bez konieczności ręcznego pobierania zbioru danych w celu zapewnienia powtarzalności eksperymentów.
GET(url, write_disk(data_file <- tempfile(fileext = ".xlsx")))
raw_df <- read_excel(data_file)
Zbiór danych składa się z 6120 rekordów medycznych opisywanych przez 81 atrybutów, z których pierwsze 7 odnosi się bezpośrednio do samego pacjenta (jego identyfikator, data przyjęcia do szpitala, wiek, …), a pozostałe 74 atrybuty zawierają informacje o wynikach przeprowadzonych badań.
Zbiór danych zawiera informacje o 375 pacjentach, u których stwierdzono obecność choroby COVID-19. Dane zostały zebrane na przestrzeni od 2020-01-10 do 2020-03-04.
Wśród pacjentów najliczniejszą grupę pod względem płci stanowią mezczyzni - 224 pacjentów (59.7%).
Średnia wieku pacjentów wynosi 58.83, pełne informacje o rozkładzie wieku pacjetnów dostępne są w tabeli poniżej.| Min | Median | Mean | Max |
|---|---|---|---|
| 18 | 62 | 58.83 \(\pm\) 16.46 | 95 |
Dodatkowo do danych o pacjentach dodany został nowy atrybut wyliczeniowy - age_group przypisujący pacjentów do dziesięcioletnich grup wiekowych, powstały na podstawie podziału atrybutu numerycznego - age. Najliczniejszą grupą wiekową wśród pacjentów jest (60,70] - 114 pacjentów (30.4%), natomiast najmniej liczną grupą jest (10,20] - 2 pacjentów (0.5%).
Statystyki dotyczące demografii pacjentów zostały również zwizualizowane na wykresach kołowych poniżej:
Minimalnym czasem hospitalizacji był 0:02:01:57 (dd:hh:mm:ss), a maksymalnym - 35:04:05:54. Średnim czasem jaki pacjenci spędzili w szpitalu był 10:20:29:08.
Pełna dystrybucja czasu hospitalizacji jest przedstawiona w tabeli poniżej:
| age_group | Mean | Max | Min | Standard_deviation |
|---|---|---|---|---|
| (10,20] | 8:16:05:49 | 14:15:13:49 | 2:16:57:48 | 8:10:25:16 |
| (20,30] | 11:13:04:00 | 17:16:22:40 | 0:10:44:51 | 4:22:31:21 |
| (30,40] | 14:02:16:38 | 35:04:05:54 | 0:03:50:46 | 7:02:59:42 |
| (40,50] | 13:13:26:15 | 29:04:45:21 | 0:12:57:09 | 7:12:09:37 |
| (50,60] | 11:07:50:50 | 33:17:10:13 | 0:03:40:29 | 7:16:55:57 |
| (60,70] | 10:10:54:04 | 29:00:38:10 | 0:03:55:53 | 7:13:14:12 |
| (70,80] | 7:18:15:24 | 27:18:53:29 | 0:11:46:57 | 6:14:10:13 |
| (80,90] | 9:00:41:43 | 31:16:54:28 | 0:02:01:57 | 8:04:20:22 |
| (90,100] | 4:10:55:24 | 4:18:31:13 | 3:21:22:44 | 0:09:38:49 |
Pacjenci byli poddawani badaniom ze średnią częstotliwością 2.26 badania na dzień. Największą ilością przeprowadzonych badań podczas pobytu pacjenta w szpitalu było - 59, natomiast najmniejszą - 1.
Bezwględna śmiertelność pacjentów ze zbioru danych wyniosła 46.0%.
Tabela z procentowo wyrażoną śmiertelnością w zależności od grupy wiekowej| survived | (10,20] | (20,30] | (30,40] | (40,50] | (50,60] | (60,70] | (70,80] | (80,90] | (90,100] | Combined |
|---|---|---|---|---|---|---|---|---|---|---|
| yes | 50.00% | 100.00% | 97.73% | 81.08% | 65.15% | 44.44% | 12.96% | 14.81% | 0.00% | 54.02% |
| no | 50.00% | 0.00% | 2.27% | 18.92% | 34.85% | 55.56% | 87.04% | 85.19% | 100.00% | 45.98% |
Najbardziej narażeni byli pacjenci z grupy wiekowej (90,100], których śmiertelność wyniosła 100.0%, ale jest to również jedna z najmniej licznych grup pacjentów.
Najmniej narażoną grupą byli pacjenci w wieku (20,30], u których nie zanotowano przypadków śmiertelnych.
Dodatkowo ze względu na małą reprezentację grupy wiekowej (10,20], tylko dwóch pacjentów - przypadek śmiertelny oraz ozdrowiciel, współczynnik śmiertelności (oraz przeżywalności) jest przekłamany.
Pozostałe atrybuty zbioru danych (znajdujące się w kolumnach od 8 do 81) zawierają wyniki badań pacjentów. Podstawowe statystyki dotyczące biomarkerów znajdują się w tabeli poniżej:
Najczęściej mierzonymi czynnikami podczas przeprowadzonych badań były: liczba leukocytów (White blood cell count) oraz liczba erytrocytów (Red blood cell count).
| Czynnik | Liczba badan | % badan |
|---|---|---|
| Interleukin 10 | 267 | 4.36% |
| Interleukin 2 receptor | 268 | 4.38% |
| Interleukin 8 | 268 | 4.38% |
| Tumor necrosis factora | 268 | 4.38% |
| Interleukin 1ß | 268 | 4.38% |
| Interleukin 6 | 272 | 4.44% |
| Czynnik | Liczba badan | % badan |
|---|---|---|
| Red blood cell count | 1127 | 18.42% |
| White blood cell count | 1127 | 18.42% |
| Serum potassium | 980 | 16.01% |
| calcium | 979 | 16.00% |
| serum sodium | 975 | 15.93% |
| Serum chloride | 975 | 15.93% |
Poniższa macierz zawiera wartości korelacji pomiędzy atrybutami biomedycznymi testowanych wśród pacjentów. Istnieje kilka korelacji pomiędzy atrybutami których wartość wynosi 1. Dzieje się tak ze względu na dwa czynniki: duża liczbę przeprowadzonych badań i mała oraz zmienna liczba sprawdzanych czynników podczas pojedynczego badania. Przykładowo korelacja czynników hemoglobin oraz HCV antibody quantification wynosi 1, ale liczba rekordów w jakich te dwa czynniki były wspólnie mierzone wynosi 2, przez co wartość ich korelacji może być błędnie zawyżona.
W artykule powiązanym z analizowanym zbiorem danych - Yan, L., Zhang, HT., Goncalves, J. et al. An interpretable mortality prediction model for COVID-19 patients - przedstawionych zostało 10 najważniejszych atrybutów zbioru danych na podstawie wyników modelu klasyfikatora Multi-tree XGBoost. Trzema najważniejszymi wyselekcjonowanymi cechami były (w kolejności malejącej ważności atrybutów): Lactate dehydrogenase (dehydrogenaza mleczanowa, LDH), lymphocyte count (liczba limfocytów) i High sensitivity C-reactive protein (wysoko czułe białko C-reaktywne, hs-CRP).
Poniżej przedstawiony został wykres, pochodzący z wymienionego artykułu, przedstawiający 10 najbardziej kluczowych atrybutów zbioru danych przy szacowaniu śmierelności pacjentów.
Autorzy m.in. poniższych publikacji:
potwierdzają obserwacje, iż badanie poziomów: enzymu dehydrogenazy mleczanowej (LDH), wysokoczułego białka C-reaktywnego (hs-CRP) oraz liczby limfocytów pacjenta mogą służyć przy szacowaniu śmierelności/dotkliwości pojedynczych przypadków koronawirusa.
Na podstawie opisanej wyżej oceny ważności atrybutów, wybranymi atrybutami do przygotowania modelu klasyfikacji zostały: Lactate dehydrogenase (dehydrogenaza mleczanowa, LDH), lymphocyte count (liczba limfocytów) i High sensitivity C-reactive protein (wysoko czułe białko C-reaktywne, hs-CRP).
W zbiorze danych występują uszkodzone rekordy - są to wiersze, w których nie ma żadnych wyników badań pacjentów oraz brakuje daty wykonania badania (kolumna test_date, RE_DATE w nieprzetworzonym zbiorze danych). Rekordy uszkodzone zostały usunięte w przetwożonym zbiorze danych.
Liczba rekordów przed usunięciem - 6120, liczba rekordów po usunięciu - 6106.
Ze względu na dużą liczbę badań oraz na to, że nie każde badanie zawierało taki sam zestaw mierzonych czynników zbiór danych został pogrupowany ze względu na pacjentów, a wartości parametrów: LDH, limfocyty i hs-CRP zostały wypełnione średnimi zmierzonymi wartościami niniejszych atrybutów dla poszczególnych pacjentów.
Zbiór danych został podzielony w proporcji 75% danych przydzielonych do zbioru uczącego (parametr p), natomiast podanym atrybutem do stratyfikacji (parametr y) był atrybut decyzyjny survived.
inTraining <- createDataPartition(
y = patients_results$survived,
p = .75,
list = FALSE)
Dane zostały poddane powtarzanej kroswalidacji oraz podany został zakres parametru mtry którego wartości będą służyły do optymalizowania klasyfikatora.
rfGrid <- expand.grid(mtry = 10:30)
gridCtrl <- trainControl(
method = "repeatedcv",
summaryFunction = twoClassSummary,
classProbs = TRUE,
number = 10,
repeats = 5)
Wybranym modelem klasyfikatora wykorzystywanym do predykcji śmiertelności pacjentów został Random Forest. Wartość parametru mtry=29 została dobrana automatycznie podczas procesu uczenia.
set.seed(23)
fitTune <- train(survived ~ .,
data = training,
method = "rf",
metric = "ROC",
preProc = c("center", "scale"),
trControl = gridCtrl,
tuneGrid = rfGrid,
ntree = 10)
fitTune
## Random Forest
##
## 264 samples
## 4 predictor
## 2 classes: 'no', 'yes'
##
## Pre-processing: centered (4), scaled (4)
## Resampling: Cross-Validated (10 fold, repeated 5 times)
## Summary of sample sizes: 237, 237, 238, 238, 238, 238, ...
## Resampling results across tuning parameters:
##
## mtry ROC Sens Spec
## 10 0.9759524 0.9366667 0.9456190
## 11 0.9744603 0.9333333 0.9512381
## 12 0.9735992 0.9333333 0.9500952
## 13 0.9753849 0.9383333 0.9484762
## 14 0.9757579 0.9333333 0.9458095
## 15 0.9756548 0.9350000 0.9540000
## 16 0.9745833 0.9400000 0.9459048
## 17 0.9747381 0.9400000 0.9457143
## 18 0.9745397 0.9400000 0.9470476
## 19 0.9760992 0.9400000 0.9454286
## 20 0.9762698 0.9383333 0.9470476
## 21 0.9762976 0.9316667 0.9499048
## 22 0.9731310 0.9266667 0.9441905
## 23 0.9742063 0.9366667 0.9526667
## 24 0.9730635 0.9316667 0.9482857
## 25 0.9746667 0.9333333 0.9498095
## 26 0.9773135 0.9366667 0.9526667
## 27 0.9727619 0.9383333 0.9495238
## 28 0.9765516 0.9433333 0.9513333
## 29 0.9766071 0.9383333 0.9512381
## 30 0.9741667 0.9316667 0.9484762
##
## ROC was used to select the optimal model using the largest value.
## The final value used for the model was mtry = 26.
Poniżej przedstawiona została macierz pomyłek wygenerowana dla wyuczonego klasyfikatora. Miara false positive dla tego klasyfikatora wynosi - 0.05, a miara false negative wynosi - 0.9787234
| no | yes | |
|---|---|---|
| no | 38 | 2 |
| yes | 1 | 46 |
Dokładność klasyfikatora wynosi - 96.55%, a natomiast miara Kappa przyjmuje wartość - 93.05%.
| score | |
|---|---|
| Sensitivity | 0.9743590 |
| Specificity | 0.9583333 |
| Pos Pred Value | 0.9500000 |
| Neg Pred Value | 0.9787234 |
| Precision | 0.9500000 |
| Recall | 0.9743590 |
| F1 | 0.9620253 |
| Prevalence | 0.4482759 |
| Detection Rate | 0.4367816 |
| Detection Prevalence | 0.4597701 |
| Balanced Accuracy | 0.9663462 |